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Abstract. We present a set of well-posed constraint-preserving boundary conditions for a 
first-order in time, second-order in space, harmonic formulation of the Einstein equations. The 
boundary conditions are tested using robust stability, linear and nonlinear waves, and are found 
to be both less reflective and constraint preserving than standard Sommerfeld-type boundary 
conditions. 
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1. Introduction 

The numerical solution of the Einstein equations has, especially in the last two years, become 
a successful means of tackling many significant physical questions. The most topical of these 
questions concern the simulation of potential sources for gravitational-wave detection, and the 
propagation of gravitational waves. These problems are most commonly solved using finite- 
size computational domains, and this involves imposing a boundary condition on the physical 
system being simulated. 

The standard treatment is to place a time-like boundary at a fixed coordinate location, 
and impose boundary conditions on the dynamical variables there. The particular conditions 
that are enforced ideally satisfy a number of properties. Most importantly, in order to 
ensure stability of the system, they should be compatible with the interior evolution equations 
so that the discretised system forms a well-posed initial-boundary-value problem (IBVP). 
Secondly, they should take into account the fact that Einstein evolutions always involve 
constraint equations as well as time evolution equations, and satisfy the constraints at all times. 
Otherwise, constraint violations introduced by the boundaries are likely to drive the evolution 
away from an Einstein solution. Finally, the boundary conditions should be compatible with 
physical considerations affecting the accuracy of the solution: they should be transparent to 
outgoing radiation, and restrict the amount of spurious incoming radiation from beyond the 
computational domain, which is assumed to contain all of the dynamics of interest. 

To date, only the system of Friedrich-Nagy fl] satisfies the above conditions for the 
fully nonlinear vacuum Einstein equations. In this formulation the evolution equations 
are expressed first-order symmetric hyperbolic form and maximally dissipative boundary 
conditions guarantee well-posedness. There have been several attempts to approach the initial 
boundary value problem for the linearized Einstein equations. More recently, Kreiss and 
Winicour proposed well-posed, constraint-preserving boundary conditions for the linearized 
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Einstein equations using the principle of frozen coefficients and pseudo-differential theory of 
systems for the first order system, which they then extrapolate to second-order 0. Other 
approaches have been suggested by Rinne for non-reflecting boundary conditions which 
control incoming radiation by specifying data for the incoming fields of the Weyl tensor J3), 
and by Buchman and Sarbach, who have followed a similar route but specifying the incoming 
fields at the boundary J4) . 

The approach which we introduce in this paper is partially derived from a method 
first discussed in a series of related papers by Kreiss, Winicour and collaborators |2][5][6], 
combined with the SBP energy method discussed in Refs. J7] [8] [9). By deriving energy 
estimates for the semi-discrete system using the "summation by parts" rule (defined below), 
one can ensure well-posedness |[T0l [TT1 [8] [T2]. By applying this approach to boundary 
conditions which are radiation controlling and constraint-preserving, we are able to construct 
an IB VP which satisfies all of the above conditions in the linearized regime. 

The conditions are derived for a harmonic formulation of the Einstein equations which 
has been implemented in lfT3l[T4l . The evolution equations of the formulation, given explicitly 
in the next section, are first-order in time, second-order in space. We approximate these 
equations using standard finite-difference techniques, however to ensure a well-posed discrete 
IBVP, we have worked out finite-difference operators for this system which satisfy the 
summation by parts property. Since our computational domain uses Cartesian coordinates 
on a cube, we have had to develop consistent operators for the corners and edges, as well. 
Following the developments of ||2] Q3] and |[T6l [TTl . we are able to construct boundary 
conditions of a Sommerfeld type, which are both well-posed and satisfy both the Einstein 
and Harmonic constraints. 

We have used the newly constructed boundary conditions in a number of practical tests 
and found them to perform extremely well in comparison with other standard techniques. Test 
evolutions include linear and nonlinear waves. In each case, the new boundary conditions are 
found to be more transparent to outgoing waves, as well as reducing the overall constraint 
violations on the grid. Further, the evolutions are stable against perturbations by high- 
frequency constraint violation ("noise") added to the data, providing a strong demonstration of 
their robustness. Tests were also done for black hole space-times. For head-on collisions and 
inspiral, our boundary conditions showed improvements in reducing reflections and constraint 
preservation, and thus improved the waveform accuracy, but as the standard treatment was not 
long term stable for these tests due to instabilities at the excision boundary, we did not feel 
that it was appropriate to display in our results. 

The plan of the paper is as follows: in the following section we introduce our harmonic 
evolution system, and describe the evolution variables, equations, and constraints. Section|2] 
describes the implementation of a harmonic evolution system for the Einstein equations. In 
particular, Section [2721 describes the construction of finite-difference operators which ensure 
that the discretisation of our evolution equations remains well-posed, including the boundary 
faces, corners and edges. The boundary treatment is described in Section[5] In Section llOl we 
present the derivation of constraint preserving Sommerfeld type conditions for this system. 
Finally, in Section [4] we discuss a number of test cases to which these boundary conditions 
have been applied, and demonstrate their usefulness in ensuring stability and improving 
accuracy in a variety of scenarios. 
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2. The harmonic evolution system 

2.1. Formulation of the evolution equations 

The decomposition of the Einstein tensor into evolution equations and constraints leaves 
four degrees of freedom in the space-time metric that are not set by the field equations 
themselves, but can be freely specified. In a 3 + 1 approach, these four degrees of freedom are 
determined by the choice of the lapse and shift, which amounts to specifying four out of ten 
metric components. The Arnowitt-Deser-Misner ("ADM") equations lfl8l are a well known 
reduction of the Einstein system corresponding to this style of gauge choice. 

An alternate approach to fixing the gauge degrees of freedom specifies the action of 
the wave operator on the coordinates, regarded as four scalar quantities. This is done by 
first choosing four functions F a and then constructing a coordinate map x a subject to the 
condition lfl9l that the d' Alembertian of each coordinate is 

Dx a = -^=3^ (V^99^d x a ) = F a , (1) 

rewriting Eqs. ([TJ as constrained variables 

C a := Ux a - F a = , (2) 

and using them in combination with the Einstein tensor one obtains the generalized 
harmonic evolution system 

:= G"" - X7^C U) + ig Ml/ V Q C* Q = . (3) 

In terms of these variables, the vacuum Einstein equations are a system of ten wave equations 
acting on the metric components, coupled through the coefficients of the wave operator and 
the source terms. 

Using the densitized inverse metric g^ v := \f—gg llv as evolution variables, the harmonic 
constraints (f2]i take the form 

C a = ^d8g ap — F a — , (4) 

while for the evolution equations we obtain 



-g 

+±gi»( 1 J!? = (d p g)(d t ,g) + V~gT T pa d T g^ + -^{d a g)d p g^ 

+2V^gV (AI F^ - y/^gg^VpF" + ^J~gA^ = , (5) 

where in the final term we have allowed for a constraint adjustment function which may 
depend on the metric and its first derivatives, 

A^:=C^(/, Vi 9 tV ). (6) 

The constraint adjustment implemented in the code are given from l20l and have the form 

All « _j* CPdr » + fl2CP y c»C v - -^=C^k, (7) 

^ py £ + e aT C°C T ^Zgli 

where the a; > are adjustable parameters, e aT is the natural metric associated with the 
Cauchy slicing, and e is a small positive number chosen to ensure regularity. 
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Assuming that the gauge source functions F a are also chosen such that they do not 
depend on derivatives of the metric, then the principle part of Eq. (0) consists of only its first 
term. That is, we have a set of ten wave equations of the form 

d P (9 pa d a r v ) = , (8) 

where S^ v are non-principle source terms consisting of at most first derivatives of the 
evolution variables. By implication, this system inherits the property of the well-posedness of 
the initial-boundary value problem for the wave equation. 

It is essential to have all of the initial data constructed in a way that satisfies the conditions 

C = 0, d t C p = 0, (9) 

as well as a construction of the boundary data that implies a homogeneous boundary condition 
for the constraints. However, by satisfying these conditions, we arrive at a well-posed 
IBVP for the constraint propagation system. Recent work by Kreiss, Winicour, Reula and 
Sarbach ifTTl demonstrates that it is possible to construct such boundary data while keeping the 
IBVP of the evolution system of the metric variables well-posed. In the following sections we 
implement and test such boundary conditions and compare them with simpler (unconstrained 
SAT and non-SAT) boundary treatments for a number of test-problems. 

In order to understand the feasibility of Eq. §5$ as an unconstrained evolution system, 
one needs to have insight into the associated constraint propagation system 



DC" = S"{g, dg, d*g, C, dC, A, OA) , (10) 

where S p is a source term dependent on the metric, the constraints, the constraint adjustment 
term, and their derivatives. 

The principal part of Eq. (ITOb is, again, that of a wave operator, implying the connection 
to results regarding the well-posedness of the initial-boundary value problem of the wave 
equation. 

We have implemented the generalized harmonic evolution system (0, cast in a form 
that is first differential-order in time, and second-differential order in space. The auxiliary 
variables 

QT = nPd p ~g a0 , (11) 

are used to eliminate the second time-derivatives, where n p is time-like and tangential to the 
outer boundary lfl3l . The resulting evolution system takes the form 

fitsT = -^sT + ^Q^, (12) 

9 9 

dtQT = -d t ((gv - dfg^ - di (CQ^ + S^(9, dg, F, OF) , (13) 

where S^ v {g, dg, F, dF) are non-principle source terms consisting of at most first derivatives 
of the evolution variables and are determined by our choice of gauge. 



2.2. Discretisation and finite differencing 

The numerical implementation of (1121113b follows the "method of lines" approach, which 
applies to systems which can be cast in the form of an ordinary differential equation containing 
some spatial differential operator L 

d t q = L(q). (14) 
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The time integration can be carried out using standard methods, such as the Runge-Kutta 
algorithm. 

Spatial derivatives on the right-hand-sides of (11211131 ) are computed by finite differencing 
on a uniformly spaced Cartesian grid. We have implemented finite difference stencils which 
are fourth-order accurate over the interior grid and second-order accurate at the boundaries. 
To ensure well-posedness of the semi-discrete system, we need to obtain an estimate on the 
energy growth of the system. To do this, we have used difference operators D which satisfy 
the "summation by parts" (SBP) property. A discrete operator is said to satisfy SBP for a 
scalar product E = (u,v) — J^u ■ vdx if 

(u,Dv) + (v,Du) = (u-v)\ b a , (15) 

holds for all functions u, v in the domain [a, b] . This is the discrete analog of the integration by 
parts property of continuous functions. By integrating for our energy estimate using the SBP 
property of our difference operators, we ensure that boundedness properties of the continuum 
energy estimate carry over to the discretised system. We can construct these difference 
operators, including numerical boundary conditions in a consistent way, for the system of 
equations in dH). 

We follow the procedure outlined by Strand |iTT]| in constructing finite difference stencils 
D of a given order, r, such that 

nil 

Du = — + 0{h T ) 1 (16) 
dx 

and which satisfy the SBP property (TT3T >. Briefly, given a state vector u = (uq, u\, . . . , u n ) T 
on ?i grid points, we construct a finite difference operator D as a matrix acting on u. The 
coefficients of D can be represented as products of the standard operators 

Dqx fi.j.k — ~^2}l ^f^^'i'^ Ijjjfc) ' 

D+x fij.k ~T {fi+l,j,k fij,k) j 

D-xfi,j,k = T (fi,j,k — fi-l,j,k) ■ (17) 

They are determined up to the boundaries of the domain by solving the set of polynomials 

dr m 

Dx m - — — =0, m = 0, 1, . . . , r, (18) 

dx 

which establish the order of accuracy r of the approximation. The SBP rule ( TT~5T > provides an 
additional set of restrictions, 

(u,Du) =-iu 2 (0), (19) 

and 

(u + v,D(u + v)) h = (D (u + v) ,u + v)h - (u + v ) 2 , (20) 

which should hold for all u, v in the half line divided into intervals of length h > 0. Following 
Strand ifTTI . we can solve these conditions explicitly for the stencil coefficients of the first 
derivative operator D. It is trivial to obtain a second derivative operator simply by repeated 
application of the derived first derivative operator. However, this results in a very wide and 
thus impractical stencil, and instead we use the second derivative SBP operators described 
in 11241 [121 . The explicit expressions for the finite difference stencils which we use are given 
inlH. 
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The above considerations apply to the construction of difference operators along a single 
coordinate direction. We can derive a 3D SBP operator by applying the ID operator along 
each coordinate direction. It can be shown that the resulting operator also satisfies SBP with 
respect to a diagonal scalar product 

(u,v) H = h x h y h z ^a i (7 j (j k u ljk ■ v ljk , (21) 

ijk 

where <7j, a k are the coefficients of the corresponding inner product in each of the 
coordinate directions. The norm H is defined such that for a discrete inner product (u, v)h = 
u T Hv, where H — H T > 0. Note that this is only true if the norm, H, is diagonal. Here we 
restrict ourselves to this case. 



3. Boundary Treatment 

3.1. Well-posed Boundary Conditions 

We have constructed finite differencing operators which satisfy summation by parts, and thus 
can use the rule dT~5b as a tool for deriving an energy estimate and ensuring well-posedness 
of the semi-discrete system. For the continuum system, we have a well defined energy 
estimate which can be used to bound solutions. Through use of the SBP-compatible derivative 
operators defined in the previous section, we ensure that an energy estimate also holds for the 
semi-discrete system. If this energy estimate bounds the norm of the solution in a resolution 
independent way, then we have a stable semi-discrete system. Optimally, we would like the 
norm of the semi-discrete solution to satisfy the same estimate as the continuum solution. 
To establish well-posedness we impose boundary conditions based upon the energy norm 



£ = \\u{t, .)\\ z = (u,u) = / u-Hudx (22) 

Jn 

where u(t, .) is the solution of the IB VP at time t, and H is a symmetric positive definite 
matrix on the bounded domain fi. We require that 

£ {t) < C{t)£ (0) , t>0, (23) 

with C (t) independent of the initial and boundary data, so that the solution is bounded by the 
energy at time t = for all t. 

As an instructive example which contains the essential features of the derivation for the 
Einstein equations, we derive explicitly the energy estimate for the wave equation with shift 
in the Appendix. 

We require that the energy, £^ n > = \\u (•, t) || 2 , satisfies d23l for positive times, that is, for 
the duration of a simulation the energy is bounded. The use of simultaneous approximation 
terms (the SAT or 'penalty') allows us to choose values for the free parameters in the boundary 
terms which conserve the energy in the system. We determine the time dependence of the 
energy for this system in order to derive coefficients for our penalty terms at the boundary 
points which give a well-posed semi-discrete system. For the wave equation with shift the 
semi-discrete evolution equations which results from the SBP-S AT calculation in the appendix 
(Sec. [5]) are 

u tt = - ^H^D^ut - ^ff-^g'u - -^ r H~ 1 E Qz (a 0z u t + (3 0l S lU + 6 0l u) 
H E Ni (a Ni u t + flNiSiii + 5 0i u) , (24) 



7 tt 0N„ 
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where is the discrete first derivative operator in the i direction and is the discrete 
second derivative operators which are blended to sideways differencing near the boundaries. 
This equation, as a result of the application of the SAT terms, satisfies the energy conservation 
equation d£/dt = 0. The corresponding calculation for the Einstein equations, Eq. ( fT2iJT3l ) 
mirrors this calculation, except with the inclusion of source terms which do not themselves 
modify the boundary treatment. 



3.2. Constraint preservation 

In ref. (T3|, we used a somewhat ad-hoc boundary condition, which applies a Sommerfeld-like 
dissipative operator to all ten components of the metric 

3 t + dx--W"-<?n = 0- (25) 



This follows the physically motivated reasoning that far away from a source, the evolution 
variables each satisfy a generally radial outgoing wavelike behaviour. The condition is 
particularly simple to apply, and has been used extensively in evolutions using a conformal- 
traceless formulation of the Einstein equations (see, for example, 1251 ). where the choice 
of evolution variables has so far hindered the development of a more rigorous boundary 
treatment. In fact, in simulations where the boundaries have been pushed to large distances 
(for instance through the use of mesh refinement), the condition has proven to be useful 
enough to allow for long-term stable evolutions. Eventually, however, boundary effects do 
contaminate the interior grid, and can lead to a loss of convergence or the accuracy required 
to resolve delicate physical features. The conditions given by Eq. d25l l make no effort to 
satisfy the Einstein constraints, and thus can over time drive the solution away from a solution 
of the full Einstein equations. 

For the Einstein equations in harmonic form, it is possible to derive consistent boundary 
conditions by explicitly evaluating the constraint propagation system. This has been done for 
the first order harmonic evolution system described by Lindblom et al. |26l . who have derived 
consistent conditions based on limiting incoming characteristics. 

Alternatively, Kreiss and Winicour [2] have demonstrated a set of Sommerfeld type 
boundary conditions, which are strongly well posed, as well as preserving the harmonic 
constraints. The well-posedness follows from results in pseudo-differential theory of strongly 
well-posed systems, and applies to a broad class of conditions. We can apply their results 
directly to the generalized harmonic evolution system used here. The harmonic constraints, 
Eq. ©, provide conditions for the time components of the metric: 

- dtg^ - d x g" x - d v g m - d z g" z - = . (26) 

The remaining metric components are determined by applying the Sommerfeld-type 
condition, Eq. ([25), in a hierarchical fashion, using previously determined components as 
required: 

d x + d t + ^j(g AB -g AB )=0, (27) 
d x + d t + -) (g tA - g xA - g tA + g xA ) = , (28) 



T 



(d x + 9 t + ^j {g u - ig xt + g xx - fl? + 2gf - = . 



(29) 
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These particular conditions are chosen to ensure well-posedness of the solution, but are not 
unique. They lead to the following explicit conditions on the positive x boundary: 

d x g°» - x g^ - dyg*" - d z9 3 » - F» , (30) 

1 



(0 x + t )g^ 
(d x + d t )g n 



(d x + dt )(2g '- g 00 )-^(g n -2g '+g m ) 

r 



Or 



(0 x +0 t )g 1A = (0 x + t )(g 

,1A 



OA 



9° A ) 



1 



(0 x +0 t )g 



AB 



-(9 AB 



9l A ) 



~0A 



9°o A ) 



x g 



1A 



9 AB )+0 x g£ 



AD 



(31) 

(32) 
(33) 



We combine the results of the previous section (see Appendix) with these constraint 
preserving conditions, to arrive at expressions for the evolution equations for from 
Eq. ( fT3l with the new penalties derived in the appendix and shown in Eq. 
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(34) 



(35) 



where the p^ u are determined by Eqs. d30ll-(l33]l. For example 

P °x = Si+g * - (Si+g* + D A+ r A + D B+ r B + F») , 

corresponds to the constraint conditions in Eqs. (f30b , where i is the direction outward from 
the boundary face, Si± is the stencil for sideways finite differencing on the boundary, and A, 
and B are tangent to the face. 



4. Applications 

The boundary prescription described in the previous section has been implemented for our 
harmonic Einstein evolution code ( |[T3l and Sec. [2]). We have carried out tests comparing 
three boundary configurations. The first, which we refer to as "standard Sommerfeld" simply 
applies Eq.|25]to each evolution variable on each face of the cubical evolution domain, which 
was the boundary implementation used in [13]. The second ("SAT") applies the boundary 
treatment derived in Sec. 13.11 and the third ("CP-SAT") improves on this by implementing 
the constraint preserving conditions of Sec. 13.21 We find that in each case, the SAT and CP- 
SAT boundary conditions respectively improve on the standard Sommerfeld condition in their 
ability to reduce boundary reflections and constraint violations over time. 



4.1. Shifted waves 

As a first test of the methodologies outlined in the previous section, we consider a simplified 
non-relativistic example problem which demonstrates the effectiveness of the SAT method. 
One of the challenges of designing boundary treatments which control the energy growth for 
black hole space-times in commonly used gauges is the problem of non-zero shift. A useful 
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Figure 1. The evolution of <j> for fiat-space wave equations with a constant shift in the x- 
direction. As initial data we have used a spherical Gaussian pulse of amplitude 1.0 and width 
1.0, on a grid 8 (121 grid points) units in size. Thin lines are the Sommerfeld-type boundary 
conditions without the SAT terms applied, whereas thick lines use the SAT boundary treatment 
given by Eq. jA.lOt . 



problem which has been used as a toy model for the full Einstein equations is the shifted scalar 
wave equation Il27ll20l . 



with shift vector /3 l = g lt / ' g tt (see Eq. ( IA.lt ). In the appendix, we have explicitly derived the 
boundary treatment of this problem, which has been implemented in a 3D evolution code. 

In Fig. [T] we display results from evolutions of a Gaussian wave packet, for various 
constant values of the shift. The Loo-norm of the energy of the solution is plotted as a function 
of time for evolutions using standard Sommerfeld type conditions, Eq. dZSt , and compared 
with the SAT conditions derived in Sec. 13.11 As the waveform impinges on the boundary, 
there is a certain amount of unphysical reflection, but the energy is largely removed from the 
grid in steps corresponding to the crossing time, as visible in Fig. |2] The boundary reflections 
are much lower in the case of the SAT boundary conditions, and the evolution is stable even 
to superluminal, \(3' l \ > 1, shifts suggesting that our conditions are stable even for outflow 
boundaries. 

4.2. Linear waves 

As a first test of the implementation of the constraint preserving boundary conditions for the 
full Einstein equations, we have considered low amplitude wave solutions of the linearized 
Einstein system. These solutions exhibit non-trivial dynamics which exercise the boundaries, 
but for which the source terms of the Einstein equations are negligible. The particular initial 
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Figure 2. The same as in Fig.[T]but shown in a logarithmic scale for 1 1 <J>oo 1 1 oo and on a longer 
timescale. Note that standard Sommerfeld boundary conditions are unstable for |/3 l | > 1. 
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Figure 3. The L2-norm of the Hamiltonian constraint for a Teukolsky wave, comparing our 
constraint-preserving boundary conditions with the standard non-SBP Sommerfeld conditions, 
as well as the purely Sommerfeld SAT algorithm to ensure well-posedness. 
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Figure 4. The tt component of the metric for a Brill wave of amplitude a = 0.5, comparing 
constraint-preserving boundary conditions with the standard Sommerfeld conditions. The 
above plot shows a two-dimensional cut in the xy plane at various times. On the right is 
the evolution of the brill wave with constraint-preserving SAT and on the left is the same 
simulation but with standard Sommerfeld type boundary conditions. 



data which we use are the quadrupole Teukolsky waves l28l . which have been used as a 
testbed in a number of numerical studies ll29l[30l[3Tl[32ll . The particular solution which we 
use follows Eppley 13311 in combining in-going and outgoing wave packets so as to produce a 
solution which is regular everywhere in the space-time. 

The overall behaviour of the evolutions using our three boundary conditions is 
summarized in Fig. [3j which plots the evolution of the L2-norm of the Hamiltonian constraint 
as a function of coordinate time, for a wave of amplitude 0.001. In each case, there 
is a reduction of the constraint violation as the wave propagates off the grid. In the 
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standard Sommerfeld case, this quickly saturates at a level of 10~ 7 , determined by the finite 
differencing resolution. In the case of the SAT boundary conditions, however, the constraint 
violation eventually reaches machine round-off due to the constraint damping in the interior 
of the domain. This happens at a much faster rate for the explicitly constraint preserving 
condition ("CP-SAT") which introduces the modification described in Sec. 13.21 It is notable 
that in this case, the initial boundary reflection which the standard Sommerfeld condition 
shares with the simple SAT treatment, is also absent. 

4.3. Nonlinear waves 

The goal of our boundary treatment is to reduce the errors introduced into the evolution 
domain during evolutions of strong field space-times involving non-linear waves, as for 
instance, generated during binary black hole evolutions. To model this problem in a simplified 
setting which does not involve complications due to excision or interior mesh-refinement 
boundaries, we have carried out tests using the nonlinear Brill wave solutions [34). These 
solutions have been studied in a number of numerical contexts, both as testbeds, as well as 
exploring the onset of black hole formation 11351 1331 [36l [371 l38l . The initial spatial metric 
takes the form 



in cylindrical (p, cf>, z) coordinates. We choose q of the form of a Gaussian packet centered at 
the origin, 



where a is a parameter which is used to set the overall amplitude of the axisymmetric wave. 
Generally we choose a value of a = 0.5 to construct a wave which is strong, but not so as to 
evolve to a black hole. As a result, we expect the initially nonlinear solution generate waves 
which propagate off the grid leaving behind a fiat space-time. 

In Fig.|4]we show a number of frames from two evolutions, displaying the metric 7*' 
component at various time instances on a grid 7 units in size. In the right column, the standard 
Sommerfeld conditions have been used, whereas on the left we have used the constraint 
preserving SAT boundary conditions. By the second frame at t = 8, the wave pulse has 
reached the boundary, and the following frames show the reflected pulse. Qualitatively, the 
CP-SAT boundary conditions show a much smoother profile, with smaller amplitude features. 
By t = 45, the wave has left the grid in the CPSBP case, to the extent that it cannot be 
seen on the linear scale of the figure. In the standard Sommerfeld case, however, there is 
still some non-trivial dynamical evolution. A more quantitative demonstration is shown in 
Fig. [5] which plots the Zv2-norm of the harmonic constraint C° as a function of coordinate 
time for three situations: The standard Sommerfeld boundary conditions ("Sommerfeld"), the 
SAT boundary conditions developed in Sec. l3.1l ("SAT"). and the constrained version of these 
boundary conditions, following the prescription of Sec. I3.2I CCP-SAT"). In the Sommerfeld 
case, the constraint violation is entirely reflected by the grid boundaries, and the value remains 
essentially constant at its initial value throughout the evolution of the evolution, even though 
constraint damping has been used on the interior code. The SAT boundary conditions, 
however, do a much better job of removing constraint violation from the grid, showing the 
exponential decrease with time that is expected from the damped solution. The constraint 
preserving boundary conditions show the strongest damping, suggesting that the constraint 
violating modes introduced by these boundary conditions are much smaller than for the SAT 
case. The evolution of the other constraint components show the same behaviour. 
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Figure 5. The L2-norm of the harmonic constraints for a Brill wave of amplitude 
0.5, comparing constraint-preserving boundary conditions with the standard Sommerfeld 
conditions, as well as the purely Sommerfeld SAT algorithm to ensure well-posedness. 
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Figure 6. Evolution of Q 00 component of the harmonic data for a Brill wave perturbed by 
random noise of a kernel amplitude of e ± 0.075, over all the grid points. This is placed on 
top of Brill wave initial data with an amplitude of a = 0.5. 
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Figure 7. Evolution of the L2-norm of the harmonic constraints for a Brill wave (a = 0.5) 
perturbed by a checkerboard noise pattern of amplitude e ± 0.1, over all the grid points, in 
order to excite the highest frequency grid mode. 



As a final test of the stability of our boundary prescription, we have carried out evolutions 
of Brill waves for which we have attempted to excite high-frequency error modes along the 
lines of the "robust stability" test [39 , 40j . This test is a means of determining whether it is 
possible for modes of any frequency within any of the grid variables to exhibit exponential 
growth during the evolution. On a numerical grid, error modes exist at fixed frequencies, 
set by the grid resolution, and the standard test consists of perturbing each variable at each 
grid point by a small amount of randomly determined amplitude e. The effect of the random 
perturbation is to seed modes which then have the potential to grow, if the system is unstable 
at that frequency. Since being first used in ll40l and proposed as a standard testbed in ll39l . the 
test has been used in a number of applications to demonstrate well-posedness of numerical 
implementations BTll42l l3l l40ll43ll44l . In Fig.[6]we applied this test by applying some kernel 
of random data to all points including the boundary points. For the SAT methods the random 
noise gets damped and then the decay of the energy looks similar to that of the standard 
brill test in Fig. [5] For the standard Sommerfeld boundary conditions the evolution becomes 
unstable at the boundaries. 

A variant of this test recognizes that in the case of an ill-posed system, the fastest 
exponential growth will result from the highest frequency mode. On a finite-difference grid, 
the frequency of this mode is set by the grid spacing. We can excite this mode by adding 
perturbations to the data in a "checkerboard" pattern, where neighboring points receive an 
opposite perturbation of fixed amplitude e. That is, we choose 

_ J +e, for i+j + k even, 
ijk ~ \ -e, for i + j + k odd. ( ' 
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In Fig. [7] we show the evolution of the L2-norm of the C° constraint component for 
the evolution of an a = 0.5 Brill wave for which each component of the initial data has 
been modified according to Eq. d39t with e = 0.1. The two versions of the SAT boundary 
conditions prove to be rather impervious to the initial data perturbation, and display essentially 
the same behaviour as in the unperturbed case, Fig. [5] It is perhaps notable that the non- 
constraint-persevering boundary conditions show a slightly slower decay rate than for the 
non-perturbed data of Fig. [5] so that it takes more than 100 time units to reach the level of 
machine round-off, whereas the constraint preserving boundary conditions reach this level 
in essentially the same amount of time as in the unperturbed case (though with a somewhat 
different decay profile). The simple Sommerfeld boundary conditions, however, are unable to 
cope with the initial perturbation and lead to an instability on a very short timescale. 

5. Conclusions 

We have examined the initial boundary value problem for the second-order formulation of the 
Einstein equations in the generalized harmonic gauge. The system of evolution equations for 
this finite-difference harmonic code was derived in fPH where it was shown to be accurate, 
stable, and convergent for long-term evolutions of black hole space-times, such as head- 
on collisions of two black holes, isolated black holes, and binary black hole inspiral and 
merger. In this paper we described the derivation, implementation and testing of a new 
boundary treatment for this system. We demonstrated that this new treatment maintained 
the validity and convergence (to lower order) seen with the standard boundary treatments. 
We additionally show that these conditions give us greater accuracy (for all reasonable 
resolutions), improved constraint preservation, improved boundary transparency, and greater 
stability in robust stability tests. 

We implemented Sommerfeld-type boundary conditions as in Eq.d25ll. which are applied 
via the simultaneous approximation term (SAT) method to control the energy growth of the 
system, and are designed to be maximally dissipative. We then establish well-posedness 
for the semi-discrete symmetric hyperbolic evolution system via the energy method fPUl by 
bounding the energy growth of the system under the assumption that the boundaries are in the 
linearized regime. We have implemented finite-differencing stencils that obey the summation 
by parts (SBP) rule | lf\ with the diagonal norm, with minimum bandwidth second-derivative 
SBP stencils as derived in J24). These stencils give fourth-order accuracy in the interior, and 
second-order at the boundary. While the standard stencils give fourth-order everywhere, we 
show that the improved accuracy of the SBP conditions more than makes up for the loss of 
two orders of convergence. 

The stability and well-posedness of the boundary conditions has been demonstrated for a 
number of test problems: shifted scalar waves, linearized waves, nonlinear waves, and random 
and high frequency stability tests. Further improved accuracy results from incorporating the 
constraint preservation into the conditions, following the prescription of 021 16) . The boundary 
conditions are still Sommerfeld type for most metric components, but we substitute conditions 
gained from enforced preservation of the harmonic constraints. This gives us four conditions 
directly from the harmonic constraints, three from the coupling of these conditions to our 
outgoing Sommerfeld-type conditions, and the three components for the directions tangent 
to each boundary face come only from our Sommerfeld-type conditions. In Sec. |4]we show 
that, as expected, these new outgoing Sommerfeld, constraint-preserving conditions retain 
the robust stability and convergence properties of the purely Sommerfeld-SBP conditions. 
The tests also demonstrate that these new conditions lead to smaller errors in satisfying the 
constraints, and are more transparent to waves propagating through the boundaries. They 
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should thus lead to more accurate evolutions than the purely Sommerfeld SBP penalty 
boundary conditions. 

In a related study, Rinne et al. l29l have considered a number of boundary treatments 
for the case of a first-order in space harmonic formulation, including the Kreiss-Winicour J2) 
treatment adopted here for a second-order system. They find that an additional physically 
motivated condition, dt^v = 0, which aims to eliminate incoming radiation, can have 
important effects in reducing physical reflections. Similar modifications may also prove 
beneficial to the second-order system presented here, though apparent reflections from the 
outer boundary are rather small even in the case of non-linear waves studied in Sec. 14.31 

With binary black hole evolutions now extending over multiple orbits, and thus many 
crossing times on conventional computational grids, boundary effects can potentially have 
a non-trivial influence on the late-time dynamics and extracted gravitational wave signals 
from such simulations. The tests provided here, including nonlinear Brill wave evolutions, 
suggest that these methods will also be effective for isolated strong sources, and thus will also 
be appropriate for black hole space-times, though these involve a number of other technical 
considerations (such as excision) which we do not explore here. The methods can be extended 
to other formulations of the Einstein equations, provided certain hyperbolicity assumptions 
are satisfied, and we are currently pursuing improvements of other commonly used systems 
such as the conformal-traceless one employed in l25ll . 
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Appendix 

As an instructive example which contains the essential features of the derivation for the 
Einstein equations, we derive explicitly the energy estimate for the wave equation with shift, 

fl?« = - 2^W) u. (Al) 

where — is the shift f3 l , and /3' /3 J — ^ is the lapse. 

We need to ensure that the energy, £( n > = \\u (•, t) || 2 , satisfies that the energy of the 
system is bounded for the duration of the simulation. The time derivative of the energy of the 
system can be re-written in semi-discrete form as follows: 

d „ d ( „ ,,r, 

= {(u t ,u t t) + (utt,u t )) tt((ui,Uj t ) + (u it ,Uj)) . (A.2) 

In this section our notation will follow that: we will use partial derivative symbols for 
continuum equations and subscripts for semi-discrete derivatives. To ensure that this quantity 
remains bounded in the semi-discrete case, we determine the energy growth which arises 
from the application of our boundary conditions, and remove this via the simultaneous 
approximation term (SAT, or "penalty") method l24l . We use a discrete second derivative 
stencil which also obeys SBP and more accurately approximates a second derivative than the 
wide stencil created from applying our first derivative twice. 
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Since we use differencing operators which obey the SBP condition, we can make use of 
Eq. dl~5b to integrate Eq. dA.2b . For the wave equation, after some algebra, this condition 
gives 



—£=-2 
dt 



Xi=Ni 



r^tt 



Xi=Ni 



(A3) 
and 



That is, the change in energy is determined by fluxes at the boundary points, x, 

Xi N i. 

On the boundary faces, we impose a set of conditions which for the moment we write in 
a generic form 

[Px,=odt + a Xi = di + S Xi=0 ] (u - u ) =0 (A.4) 



--Nd t 



a x 



di - 8 Xi=N ] (u -u Q ) =0 



(A.5) 



in terms of parameters a, B, and 5 which are indexed according to the grid face. These are 
substituted into into the estimate, Eq. (tA.3t . leading to 



— £ = -2 
dt 



a Nl 



0. 



Xi=Ni 



00,' 



Xi=0 



(A.6) 



where if is the normal to the boundary face i, and uq are data chosen to be consistent with 
the initial data. The SAT method allows us to choose values for the free parameters in the 
boundary terms which conserve the energy in the system. We first write the original shifted 
wave equation, Eq. ( IA.1I ). in semi-discrete form, explicitly including the boundary terms: 



-^H^D^u - 2^ I H- 1 D ( i 1) u t + T 0i H- l E 0i (a 0i Ut + Qi S iU + S 0t u) 



+ r Nt H 1 E Nz (a Ni u t + I3n 1 S 1 u + S^u) 



(A.7) 



The E a are vectors of length N defined as E m = (0, ... 0, 1) T and E 0i = (1, . . . , 0) T to 
be zero everywhere except at the boundary points. S, are sideways blended finite differencing 
stencils satisfying the SBP property, as described in the previous section. 

We determine the time dependence of the energy for this new system in order to 
derive coefficients t for our penalty terms which give a well-posed semi-discrete system. 
Substituting Eq. dA.7b into Eq. dA.2b . and once again making use of the SBP property, 
Eq. < n~5b . we arrive at 



dt 



£ = 



u t E Nt u t + 2 ( T 0z a 0i + — 



u t E Qz u t 



2 T Ni Ni 77 )u t E Ni SiU + 2 t ,J3 0i + —77 )u t EotSiU 



r 



(A.8) 



The free parameters tq and tm can be used to eliminate the ujE^SiU terms, by setting 



tqPo = -7 t5 /7* 



and tn/3n 
d 



7^/7 . Then, the energy evolves according to 



dt 



£ 



2/3 



r 



8 1%t 



2_ 



NiUt 



1 rsjtt 



B-\jE 0t u t = Q. 



(A.9) 



The last equality is arrived at after some algebra, substituting the boundary conditions, 
Eq. dA.44IA.5T l, and making use of the original wave equation, Eq. (IA.U . 
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The resulting semi-discrete evolution equation is given by 

it IJ 1J 

u tt = - ^H^D^ut - 2_F-i£)g) u _ _)_H-i£ 0i ( a0iUt + faSm + 5 0i u) 
1 T T POi 

+ -f-r—H^E^iaN^Ut + p Ni S iU + <5 0l u) , (A.10) 

7 PNi 

which, as a result of the application of the SAT terms, satisfies the energy conservation 
equation d£/dt = 0. 



